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Abstract 

We revise the 6* model for the Collins-Soper-Sterman resummed form factor to 
improve description of the leading-power contribution at nearly nonperturbative 
impact parameters. This revision leads to excellent agreement of the transverse 
momentum resummation with the data in a global analysis of Drell-Yan lepton pair 
and Z boson production. The nonperturbative contributions are found to follow 
universal quasi-linear dependence on the logarithm of the heavy boson invariant 
mass, which closely agrees with an estimate from the infrared renormalon analysis. 



Transverse momentum distributions of heavy Drell-Yan lepton pairs, W, or Z 
bosons produced in hadron-hadron collisions present an interesting example 
of factorization for multi-scale observables. If the transverse momentum qx of 
the electroweak boson is much smaller than its invariant mass Q, da/dqx at 
an n-th order of perturbation theory includes large contributions of the type 
a™\n m (g^ / Q 2 ) / qj< (m = 0, 1 . . . 2n — 1), which must be summed through all 
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orders of a s to reliably predict the cross section [1]. The feasibility of all-order 
resummation is proved by a factorization theorem, first formulated for e + e~ 
hadropro duct ion [2,3], stated by Collins, Soper, and Sterman (CSS) for the 
Drell-Yan process [4], and recently proved by detailed investigation of gauge 
transformations of /c^-dependent parton densities [5,6]. 

The heavy bosons acquire non-zero qr mostly by recoiling against QCD radi- 
ation. The CSS formalism accounts for both the short- and long- wavelength 
QCD radiation by means of a Fourier-Bessel transform of a resummed form 
factor W(b) introduced in impact parameter (b) space. The perturbative con- 
tribution, characterized by b < 0.5 GeV -1 , dominates in W and Z boson 
production at all values of qr- The nonperturbative contribution from b > 
0.5 GeV^ 1 is not negligible at qx < 20 GeV in the precision measurements 
of the W boson mass M w at the Tevatron and LHC [7]. The model for the 
nonperturbative recoil is the major source of theoretical uncertainty in the 
extraction of M w from the experimental data. This uncertainty must be re- 
duced in order to measure M w with accuracy of about 30 MeV in the Tevatron 
Run-2 and 15 MeV at the LHC. The nonperturbative model presented below 
approaches the level of accuracy desired in these measurements. 

The nonperturbative component [described by the function J-Np{b,Q) given 
in Eq. (4)] can be constrained in a few experiments by exploiting process- 
independence, or universality, of J-Np{b, Q), just as the universal /^-integrated 
parton densities are constrained with the help of inclusive scattering data. The 
universality of FNp{b,Q) in unpolarized Drell-Yan-like processes and semi- 
inclusive deep-inelastic scattering (SIDIS) follows from the CSS factorization 
theorem [5]. In the study presented here, we carefully investigate agreement of 
the universality assumption with the data in a global analysis of fixed-target 
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Drell-Yan pair and Tevatron Z boson production. We revise the nonperturba- 
tive model used in the previous studies [8,9] and improve agreement with the 
data without introducing additional free parameters. Renormalization-group 
invariance requires J^Npip^Q) to depend linearly on InQ [3,4]. With our lat- 
est revisions put in place, the global qr fit clearly prefers a simple function 
^Np(b,Q) with universal InQ dependence. The new J^NpipiQ) has reduced 
dependence on the collision energy \fS comparatively to the earlier fits. The 
slope of the In Q dependence found in the new fit agrees numerically with its 
estimate made with methods of infrared renormalon analysis [10,11]. 

The function J^Npib, Q) primarily parametrizes the "power-suppressed" terms, 
i. e., terms proportional to positive powers of b. When assessed in a fit, J^Npib, Q) 
also contains admixture of the leading-power terms (logarithmic in b terms), 
which were not properly included in the approximate leading-power function 
Wip(b) [cf. Eq. (4)]. In contrast, estimates of TNp{b,Q) made in the infrared 
renormalon analysis explicitly remove all leading-power contributions from 
^Np(b,Q) [11]. While the recent studies [9,10,11,12,13] point to an approx- 
imately Gaussian form of J^NpibiQ) \FNp{b,Q) oc b 2 }, they disagree on the 
magnitude oiFxpib, Q) and its Q dependence. The source of these differences 
can be traced to the varying assumptions about the form of the leading-power 
function Wip(b) at b < 2 GeV -1 , which is correlated in the fit with JF A r P (6, Q). 
The exact behavior of W(b) at b > 2 GeV -1 is of reduced importance, as 
W(b) is strongly suppressed at such b. The new improvements described here 
(excellent agreement of J-Np{b,Q) with the data and renormalon analysis) 
result from modifications in the model for W / lp(6) at b < 2 GeV -1 . The im- 
provements are preserved under variations of the large-6 form of Wlp(6) in a 
significant range of the model parameters. 
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Our paper follows the notations in Ref. [9]. The form factor W(b) factorizes 
at all b as [2,3,4] 

W(b) = £ -L-e-sMPjix^Vsfab), (1) 

3=9,9 b 

where af^/S is a constant prefactor [4], and X\^ = e ±v Qj\fS are the Born- 
level momentum fractions, with y being the rapidity of the vector boson. The 
6-dependent parton densities Vj(x,b) and Sudakov function 



S(b,Q)= J ^ A (a s (fl))]n(^j+B(a s m 



2- 

(2) 



are universal in Drell-Yan-like processes and SIDIS [5]. When the momentum 
scales Q and b /b (where b = 2e~ 7s 1.123 is a dimensionless constant) are 
much larger than 1 GeV, W(b) reduces to its perturbative part W pert (b), i.e., 
its leading-power (logarithmic in b) part evaluated at a finite order of a s : 



<J {0) 

3=g,q 

x [C® f]Ax u b; fji F )[C® fk(x 2 ,b; » F ). (3) 



Here S P {b,Q) and [C <g> f}j (x, b; fi F ) = J2a Ix d £/£C ja {x/£, (i F b)f a {£, (i F ) are 
the finite-order approximations to the leading-power parts of S(b,Q) and 
Vj(x,b). f a (x,/i F ) is the /cy-integrated parton density, computed in our study 
by using the CTEQ6M parameterization [14]. Cj a (x, fx F b) is the Wilson coef- 
ficient function. We compute Sp(b,Q) up to 0(a 2 s ) and Cj a up to 0(a s ). 

In Z boson production, the maximum of bW(b) is located at b ps 0.25 GeV -1 , 
and W pert (b) dominates the Fourier- Bessel integral. In the examined low-Q 
region, the maximum of bW{b) is located at b ~ 1 GeV -1 , where higher-order 
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corrections in powers of a s and b must be considered. We reorganize Eq. (1) to 
separate the leading-power (LP) term Wipib), given by the model-dependent 
continuation of W pert (b) to b > 1 GeV -1 , and the nonperturbative exponent 
e -^Np(b,Q)^ absorbs the power-suppressed terms: 

W(b) = W LP {b) e~ TNP ^ Q) . (4) 

At b — > 0, the perturbative approximation for W(&) is restored: W LP — >■ W pert , 
J~np — > 0. The power-suppressed contributions are proportional to even pow- 
ers of b [10]. Detailed expressions for some power-suppressed terms are given 
in Ref. [11]. At impact parameters of order 1 GeV -1 , we keep only the first 
power-suppressed contribution proportional to b 2 : 

J- np ~ b 2 (ai + a 2 \n(Q/Qo) + a z <p(xi) + a 3 (f)(x 2 )) + ... , (5) 

where ai, a 2 , and a 3 are coefficients of magnitude less than 1 GeV 2 , and 4>(x) 
is a dimensionless function. The terms a2ln(Q/Qo) and as<f)(xj) arise from 
S(b,Q) and ln[7-j(xj, 6)] in ln[W(6)], respectively. We neglect the flavor de- 
pendence of <f>(x) in the analyzed region dominated by scattering of light u 
and d quarks. J-'np is consequently a universal function within this region. The 
dependence of J-'np on InQ follows from renormalization-group invariance of 
the soft-gluon radiation [3]. The coefficient a 2 of the InQ term has been re- 
lated to the vacuum average of the Wilson loop operator and estimated within 
lattice QCD as 0.19±g;^ GeV 2 [11]. 

The preferred J-'np is correlated in the fit with the assumed large-6 behavior of 
Wlp- We examine this correlation in a modified version of the 6* model [3,4]. 
The shape of Wlp is varied in the b* model by adjusting a single parameter 
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b m ax- Continuity of W and its derivatives, needed for the numerical stability 
of the Fourier transform, is always preserved. We set Wip{b) = W pert (b*), 
with b*(b,b max ) = b(l + 6 2 /^LJ" 1/2 - W LP (b) reduces to W pert {b) as b -> 
and asymptotically approaches W pert (b max ) as 6 — > oo. The 6* model with a 
relatively low 6 maa , = 0.5 GeV -1 was a choice of the previous q F fits [8,9]. 
However, it is natural to consider b max above 1 GeV -1 in order to avoid ad 
hoc modifications of W pert (b) in the region where perturbation theory is still 
applicable. To implement W pert (K) for b max > 1 GeV" 1 , we must choose the 
factorization scale ji F such that it stays, at any b and b max , above the initial 
scale Q ini = 1.3 GeV of the DGLAP evolution for the CTEQ6 PDF's f a (x, fi F ). 
We keep the usual choice \if = C3/b*(b,b max ), where C3 ~ bo, for b max < 
bo/ Qini ~ 0.86 GeV -1 . Such choice is not acceptable at b max > b /Qi n i, as it 
would allow < Qini- Instead, we choose \if = Cs/b*(b,bo/Qini) for b max > 
b / Qini, i-e., we substitute b /Q in i for b max in fi F to satisfy fx F > Q ini at any 
b. Aside from f a (x,/iF), all terms in W pert (b) are known, at least formally, as 
explicit functions of a s (l/b) at all b < 1/Aqcd- We show in Ref. [15] that 
this prescription preserves correct resummation of the large logarithms and is 
numerically stable up to 

bmax ~ 3 GeV 

We perform a series of fits for several choices of b max by closely following the 
previous global qr analysis [9]. We consider a total of 98 data points from 
production of Drell-Yan pairs in E288, E605, and R209 fixed-target experi- 
ments, as well as from observation of Z bosons with q T < 10 GeV by CDF 
and D0 detectors in the Tevatron Run-1. See Ref. [9] for the experimental 
references. Overall normalizations for the experimental cross sections are var- 
ied as free parameters. Our best-fit normalizations agree with the published 
values within the systematical errors provided by the experiments, with the 
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Q [GeV] 

Fig. 1. The best-fit values of a(Q) obtained in independent scans of x 2 f° r the 
contributing experiments. The vertical error bars correspond to the increase of x 2 
by unity above its minimum in each Q bin. The slope of the line is equal to the 
central- value prediction from the renormalon analysis [11]. 

exception of the CDF Run-1 normalization (rescaled down by 7%). 

To test the universality of J-'np, we individually examine each bin of Q. We 
choose J-'np = a{Q)b 2 and independently fit it to each of the 5 experimental 
data sets to determine the most plausible normalization in each experiment. 
We then set the normalizations equal to their best-fit values and examine 
X 2 at each Q as a function of a(Q). For b max = 1 — 2 GeV -1 , the best-fit 
values of a(Q) follow a nearly linear dependence on InQ [cf. Fig. 1]. The 
slope a,2 = da(Q) / d(\n Q) is close to the renormalon analysis expectation 
of 0.19 GeV 2 [11]. The agreement with the universal linear \nQ dependence 
worsens if b max is chosen outside the region 1-2 GeV -1 . Since the best-fit a(Q) 
are found independently in each Q bin, we conclude that the data support the 
universality of J-'np, when b max lies in the range 1 — 2 GeV -1 . In another test, 
we find that each experimental data set individually prefers a nearly quadratic 
dependence on b, Tnp = a(Q)& 2-/3 , with |/3| < 0.5 in all five experiments. 
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Fig. 2. The best-fit x 2 an d coefficients ai, 02, and 03 in TNp(b,Q) for different 
values of b max , C3 = 60 (stars) and C3 = 2bo (squares). The size of the symbols 
approximately corresponds to la errors for the shown parameters. 

To further explore the issue, we simultaneously fit our model to all the data. 
We parametrize a(Q) as a(Q) = ai + a2in [Q/(3.2 GeV)]+a3ln [IOOX1X2] . This 
parametrization coincides with the BLNY form [9], if the parameters are re- 
named as {gi, g 2 , (7i(?3}(BLNY) — > {a 1; a 2 , a 3 }(here). It agrees with the generic 
form of ^NpibyQ) in Eq. (5), if one identifies (j>{x) = ln(x/0.1). We carry out 
two sequences of fits for C3 = 60 and C3 = 26 to investigate the stability of 
our prescription for [i F and sensitivity to 0(a 2 s ) corrections. The dependence 
on C3 is relatively uniform across the whole range of b max , indicating that our 
choice of tip for b max > b Q /Q ini is numerically stable. 

Fig. 2 shows the dependence of the best-fit x 2 , a%, a>2, and a 3 on b max . As 
b max is increased above 0.5 GeV -1 assumed in the BLNY study, x 2 rapidly 
decreases, becomes relatively flat at b max = 1 — 2 GeV^ 1 , and grows again 
at b max > 2 GeV -1 . The global minimum of \ 2 = 125(111) is reached at 
bmax ~ 1-5 GeV -1 , where all data sets are described equally well, without 
major tensions among the five experiments. The improvement in % 2 mainly 
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Fig. 3. Differences between the measured (Data) and theoretical (Theory) cross 
sections, divided by the experimental error 5 exp in each (Q,qr) bin. The values of 
X 2 for each experiment in the two fits are listed in the legend in the same order, 
ensues from better agreement with the low-Q experiments (E288, E605, and 

R209), while the quality of all fits to the Z data is about the same. This 
is illustrated by Fig. 3, which shows the differences between the measured 
and theoretical cross sections, divided by the experimental errors 5 exp , as well 
as the values of y 2 m each experiment, in two representative fits for b max = 
0.5 GeV -1 , C3 = b (squares) and b max = 1.5 GeV -1 , C3 = 2b (triangles). 
The data are arranged in bins of Q (shown by gray background stripes) and 
qr, with both variables increasing from left to right. For b max = 1.5 GeV -1 , 
the (Data — Theory) differences are reduced on average in the entire low-Q 
sample, resulting in lower x 2 m three low-Q experiments. Two outlier points 
in the E605 sample (the first point in the second Q bin and fifth point in the 
fifth Q bin) disagree with the other E288 and E605 data in the same Q and x 
region and contribute 15 — 25 units to x 2 a t any b max . If the two outliers were 
removed, one would find x 2 /d-o.f. ~ 1 at the global minimum. 

The magnitudes of ai, a 2 , and a 3 are reduced when b max increases from 0.5 
to 1.5 GeV -1 . In the whole range 1 < b max < 2 GeV -1 , a 2 agrees with the 
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p + p -> Z° + X; Vs=1 .96 TeV; Q=M Z ; y=0 



p + Cu -> n*ir + X; Vs=38.8 GeV; Q=1 1 GeV; y=0 
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Fig. 4. The best-fit form factors bW(b) in (a) Tevatron Run-2 Z boson production; 
(b) E605 experiment. In the E605 case, bW(b) are divided by the best-fit normaliza- 
tions Nfit for the E605 data, and the form factor in the Qiu-Zhang parametrization 
[12] for bmax = 0.3 GeV -1 is also shown. 

renormalon analysis estimate. The coefficient 03, which parametrizes devia- 
tions from the linear InQ dependence, is considerably smaller (< 0.05) than 
both a\ and a 2 (~ 0.2). A reasonable quality of the fit is retained if 03 is set 
to zero by hand: x 2 increases by ~ 5 in such a fit above its minimum in the 
fit with a free 03. In contrast, x 2 increases by > 200 units if 03 = gig^ is set 
to zero at b max = 0.5 GeV -1 , as it was noticed in the BLNY study. 

The preference for the values of b max between 1 and 2 GeV -1 indicates, first, 
that the data do favor the extension of the b range where WLp{b) is approx- 
imated by the exact W pert (b). In Z boson production, this region extends up 
to 3 — 4 GeV -1 as a consequence of the strong suppression of the large-6 tail 
by the Sudakov exponent. The fit to the Z data is actually independent of 
bmax within the experimental uncertainties for b max > 1 GeV" 1 . The best-fit 
form factors bW(b) for b max = 0.5 and 1.5 GeV -1 in Z boson production are 
shown in Fig. 4(a). 
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In the \ow-Q Drell-Yan process, continuation of bW pert (b) far beyond b 
1 GeV -1 raises objections, since bW pert (b) has a maximum and is unstable 
with respect to higher-order corrections at b ~ 1.2 — 1.5 GeV -1 . The dubious 
large contributions to W pert {b) in this b region would deteriorate the quality 
of the fit. The b* prescription with b max < 2 GeV -1 reduces the impact of the 
dubious terms on W(b): for b max small enough, the maximum of W pert (K) is 
only reached at b 3> 1.2 GeV -1 , where it is suppressed by e~ :pNP ^ h,Ci \ The best- 
fit form factors for the E605 kinematics, divided by the best-fit normalizations 
of the E605 data Nf it , are shown in Fig. 4(b). 

If a very large b max comparable to 1/Kq C d is taken, Wipip) essentially co- 
incides with W pert (b), extrapolated to large b by using the known, although 
not always reliable, dependence of W pert (b) on In b. Similar, but not identical, 
extrapolations of W pert {b) to large b are realized in the models [12,13], which 
describe the Z data well, in accord with our own findings. In Z boson produc- 
tion, our best-fit a(M z ) = 0.85 ± 0.10 GeV 2 agrees with 0.8 GeV 2 found in 
the extrapolation-based models, and it is about a third of 2.7 GeV 2 predicted 
by the BLNY parametrization. Our results support the conjecture in [12] that 
a 3 is small if the exact form of W pert {b) is maximally preserved. To describe 
the \ow-Q data, the model [12] allowed a large discontinuity in the first deriva- 
tive of W(b) at b equal to the separation parameter b^ x = 0.3 — 0.5 GeV -1 , 
where switching from the exact W pert (b) to its extrapolated form occurs [cf. 
Fig. 4(b)]. In the revised 6* model, such discontinuity does not happen, and 
WLp(b) is closer to the exact W pert (b) in a wider b range at low Q than in the 
model [12]. The two models differ substantially at b w 1 GeV -1 , as seen in 
Fig. 4(b). 
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To summarize, the extrapolation of W pert (b) to b > 1.5 GeV" 1 is disfavored by 
the \ow-Q data sets, if a purely Gaussian form of J~np is assumed. The Gaus- 
sian approximation is adequate, on the other hand, in the 6* model with b max 
in the range 1 — 2 GeV -1 . Here variations in b max are compensated well by ad- 
justments in di, a 2 , and a 3 , and the full form factor bW{b) stays approximately 
independent of b max . The best-fit parameters in Tnp are quoted for b max = 
1.5 GeV" 1 as {a u a 2 ,a 3 } = {0.201 ±0.011,0.184±0.018,-0.026±0.007} GeV 2 
for C 3 = bo, and {0.247 ± 0.016, 0.158 ± 0.023, -0.049 ± 0.012} GeV 2 for 
C 3 = 2b . In Ref. [15], the experimental errors are propagated into various 
theory predictions with the help of the Lagrange multiplier and Hessian ma- 
trix methods, discussed, e.g., in Ref. [14]. We find that the global fit places 
stricter constraints on F NP at Q — M z than the Tevatron Run-1 Z data alone. 
Theoretical uncertainties from a variety of sources are harder to quantify, and 
they may be substantial in the \ow-Q Drell-Yan process. In particular, x 2 
for the \ow-Q data improves by 14 units when the scale parameter C 3 in \i P 
is increased from b to 26 , reducing the size of the finite-order W pert (b) at 
low Q. The best-fit normalizations Nf it also vary with C 3 . The dependence 
of the quality of the fit on the arbitrary factorization scale \i P indicates im- 
portance of 0(a 2 s ) corrections at low Q, but does not substantially increase 
uncertainties at the electroweak scale. Indeed, the 0(a 2 s ) corrections and scale 
dependence are smaller in W and Z production. In addition, the term a 2 InQ, 
which arises from the soft factor S(b,Q) and dominates J-'np set Q — M z , 
shows little variation with C 3 [cf. Fig. 2c] . Consequently, the revised 6* model 
with b max ~ 1.5 GeV -1 reinforces our confidence in transverse momentum 
resummation at electroweak scales by exposing the soft-gluon origin and uni- 
versality of the dominant nonperturbative terms at collider energies. 
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